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Abstract.  By  revealing  complex  fiber  structure  through  the  orientation 
distribution  function  (ODF),  q-ball  imaging  has  recently  become  a  popular 
reconstruction  technique  in  diffusion-weighted  MRI.  In  this  paper,  we  propose 
an  analytical  dimension  reduction  approach  to  ODF  maxima  extraction.  We 
show  that  by  expressing  the  ODF,  or  any  antipodally  symmetric  spherical 
function,  in  the  common  fourth  order  real  and  symmetric  spherical  harmonic 
basis,  the  maxima  of  the  two-dimensional  ODF  lie  on  an  analytically  derived 
one-dimensional  space,  from  which  we  can  detect  the  ODF  maxima.  This 
method  reduces  the  computational  complexity  of  the  maxima  detection,  without 
compromising  the  accuracy.  We  demonstrate  the  performance  of  our  technique 
on  both  artificial  and  human  brain  data. 


1  Introduction 

Diffusion-weighted  MRI  significantly  extends  the  scope  of  the  information  obtained 
from  MRI,  from  being  solely  spatially  dependent  to  being  defined  on  the  spatial- 
orientational  domain.  Fiber  microstructure  and  orientation  are  inferred  using  this 
modality  from  the  locally  measured  diffusion  profile  of  water  molecules.  Diffusion 
tensor  imaging  (DTI)  [1]  effectively  models  the  diffusion  in  single-fiber  voxels  as  a 
Gaussian  represented  by  its  covariance  tensor.  As  for  more  complex  fiber 
architecture,  q-ball  imaging  (QBI)  [2]-[6]  has  been  very  successful  in  revealing 
intravoxel  fiber  orientations  by  introducing  the  orientation  distribution  function 
(ODF)  as  the  probability  of  diffusion  in  a  given  direction. 

Contrary  to  DTI,  where  the  principal  diffusion  direction  can  be  readily  computed 
as  the  major  eigenvector  of  the  diffusion  tensor,  QBI  provides  a  continuous  spherical 
function  which,  although  clearly  illustrates  the  major  diffusion  orientations  as  its 
maxima,  does  not  directly  quantify  them.  Diffusion  directions  as  vectors  carry  less 
information  than  the  ODF  itself  does.  On  the  other  hand,  their  easy  interpretation  and 
their  application  in  tractography,  e.g.,  [7] — [9],  make  the  ODF  maxima  extraction  an 
important  post-processing  step  still  to  be  carefully  addressed.  The  number  of  peaks 
can  also  be  interpreted  as  a  measure  of  white  matter  complexity.  In  addition,  unlike 
mixture  models  that  calculate  fiber  directions  by  describing  the  diffusion  signal  as  the 
sum  of  finite  discrete  unidirectional  components,  ODF  maxima  are  computed  without 
any  assumptions  about  the  existence  of  such  components. 


Exhaustive  search  via  finite  difference  method  has  been  exploited  in  the  literature 
as  a  straightforward  approach  to  ODF  maxima  extraction  [3,10].  This  generally 
requires  a  two-dimensional  (2D)  discretization  of  the  unit  sphere,  resulting  in 
computational  complexity  that  grows  quadratically  with  the  desired  resolution. 
Numerical  optimization  approaches  such  as  gradient  ascent  [11],  Newton-Raphson 
techniques  [12],  and  Powell’s  method  [13],  have  also  been  employed.  These 
techniques  require  a  guarantee  of  convergence  and  careful  initialization  to  obtain  all 
the  maxima.  Lastly,  polynomial  based  approaches,  [14]— [16],  have  been  proposed  to 
extract  the  maxima  as  a  subset  of  the  stationary  points  of  the  ODF.  These  methods 
exploit  a  transformation  of  the  real  and  symmetric  spherical  harmonic  (RSSH)  basis 
(most  efficient  for  ODF  reconstruction  [3]),  to  the  constrained  symmetric  tensor  or 
constrained  homogenous  polynomial  bases,  resulting  in  polynomial  equations  which 
are  solved  numerically. 

In  this  paper,  we  propose  a  polynomial  based  approach  to  reduce  the  problem  of 
ODF  maxima  extraction  in  the  fourth  order  RSSH  basis,  from  a  2D  search  on  the 
sphere,  to  a  one-dimensional  (ID)  one  on  an  analytically-derived  curve.  Compared  to 
the  2D  problem,  this  approach  significantly  reduces  the  computational  complexity  of 
the  search  for  the  maxima  of  the  ODF  -  or  any  antipodally  symmetric  spherical 
function  -  without  compromising  the  precision.  Contrary  to  [14]— [16],  our  method 
works  directly  in  the  RSSH  basis  and  does  not  require  the  extra  step  of  transforming 
the  RSSH  coefficients  to  other  tensor-based  bases.  We  suggest  a  discretization 
scheme  for  the  ID  exhaustive  search,  and  show  experimental  results  on  both  artificial 
and  human  brain  data. 

We  start  Sec.  2  with  a  brief  review  of  the  RSSH  basis,  and  continue  by  describing 
our  mathematical  derivation.  Experimental  results  are  presented  in  Sec.  3. 


2  Methods 

2.1  ODF  in  Real  and  Symmetric  Spherical  Harmonic  Basis 

In  this  work,  we  use  the  estimator  derived  in  [6]  to  compute  the  ODF  in  constant  solid 
angle  (CSA).  The  original  definition  of  the  QBI  ODF  [2]  does  not  include  the 
Jacobian  factor  rz,  thereby  creating  the  need  for  normalization  and  artificial 
sharpening.  In  contrast,  the  estimator  in  [6]  is  normalized,  dimensionless,  and  has 
been  shown  to  preserve  the  natural  sharpness  of  the  ODF. 

The  spherical  harmonic  basis  is  commonly  used  for  representing  spherical 
functions  such  as  the  ODF,  allowing  for  sampling  in  any  desired  direction. 
Orthonormal  spherical  harmonic  functions  are  given  by 
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where  P{n(-')  is  the  associated  Legendre  function,  and  9  and  (p  are  standard  spherical 
coordinates.  The  assumption  of  the  ODF  being  real  and  antipodally  symmetric, 
however,  makes  the  use  of  the  RSSH  basis  [3]  more  suitable.  RSSH  functions  are 


indexed  by  a  single  parameter  j  =  1(1  +  l)/2  +  m  +  1,  corresponding  to  lj  and  mj,  as 
follows  [3]: 
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The  ODF  can  be  computed  in  this  basis  first  by  using  a  minimum  square  scheme  to 
approximate  the  signal,  and  then  by  analytically  computing  the  Funk-Radon 
transform,  [2],  following  the  method  introduced  in  [3]— [5]  for  the  original  QBI,  and 
subsequently  adapted  in  [6]  for  the  CSA-QBI. 


2.2  ODF  Maxima  Extraction 

RSSH  functions,  being  smooth,  allow  us  to  find  all  the  local  maxima  of  the  ODF 
ip(8,  (p )  as  points  satisfying  the  following  properties  (subscripts  indicate  partial 
derivatives): 


with  the  Flessian  matrix  H{d,cp') 

H(0, 4>) 
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Equations  (3)  and  (4)  guarantee  that  ( 9 ,  (p )  is  either  an  extremum  or  a  saddle  point 
of  the  ODF.  Inequalities  (5)  and  (6)  filter  out,  respectively,  the  saddle  points  and  the 
local  minima  (including  possible  negative  lobes),  leaving  us  only  with  the  local 
maxima  of  the  ODF.  The  above  expressions  can  all  be  analytically  computed  for  an 
ODF  expressed  in  the  RSSH  basis.  However,  the  main  challenge  is  to  find  the  points 
that  simultaneously  satisfy  equations  (3)  and  (4).  Once  they  are  identified,  applying 
inequalities  (5)  and  (6)  to  filter  out  undesired  points  is  trivial. 

Iterative  approaches  (e.g.,  Newton  method)  may  be  applied  to  solve  equations  (3) 
and  (4).  Yet,  being  quite  sensitive  to  the  initialization,  they  are  not  guaranteed  to 
converge  to  all  the  maximum  points.  Alternatively,  an  exhaustive  search  will  result  in 
all  the  maxima  with  an  accuracy  determined  by  the  discretization  resolution. 
Nonetheless,  with  the  ODF  being  a  2D  manifold,  the  search  space,  and  consequently 
the  computational  complexity  of  the  algorithm,  grows  quadratically  with  the  desired 
resolution. 


We  will  next  show  how  the  fourth  order  RSSH  basis  makes  it  possible  to  confine 
the  search  to  a  ID  space,  thereby  creating  an  efficient  method  to  extract  the  maxima. 


2.3  Reducing  the  Dimension  of  the  Search  Space 


Let  us  assume  that  the  ODF  has  been  approximated  in  the  fourth  order  RSSH  basis,  as 

15 
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Combining  equations  (1),  (2),  and  (8),  while  substituting  the  values  of  P;m(cos0) 
from  Table  1  leads  to 
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a5sin0),  etc.  (We  drop  the  notation 


(0)  in  the  rest  of  this  subsection.) 

We  now  attempt  to  solve  Eq.  (4)  by  deriving  Eq.  (9)  with  respect  to  (p.  We  then 
divide  it  by  sin  0  cos3  0  and  rearrange  it,  while  using  the  identity  sec2  0  =  1  + 
tan2  0,  to  obtain 


(H  +  C  —  F)^  tan3  0  +  (G  +  B  —  3E)rp  tan2  0  +  (6 F  +  C)^  tan  0 
+  (B  +  4  E)^  =  0. 


(10) 


Equation  (10)  is  a  cubic  function  of  tan  0,  and  can  be  analytically  solved,  leading 
to  a  closed-form  expression  for  0(0).'  Thus,  for  each  given  (p ,  we  obtain  one,  two,  or 
three  different  real  values  for  0  which  satisfy  Eq.  (4). 

The  curve  characterized  by  the  pair  (0 (0),  0)  (Fig.  1(b))  is  in  fact  our  new  ID 
search  space  which  contains  all  the  ODF  maxima  as  points  satisfying  equations  (3), 
(5),  and  (6)  (Fig.  l(c&d)).  The  number  of  these  maxima  does  not  need  to  be  initially 
specified,  since  it  is  automatically  determined  by  the  algorithm  and  depends  on 
various  factors,  such  as  the  number  of  real  solutions  to  Eq.  (10).  This  is  particularly 
important  in  practice,  as  different  regions  of  the  white  matter  naturally  exhibit 
different  complexity.  The  maxima  can  be  found  using  a  ID  exhaustive  search  (see 
Sec.  2.4),  which  is  considerably  faster  than  exploring  the  entire  2D  manifold  of  the 
ODF.2 


1  Each  solution  of  tan  0  corresponds  to  a  unique  value  of  0  £  [0,  n).  Please  note  that  this 
approach  can  also  be  applied  in  the  RSSH  basis  of  higher  orders,  with  the  difference  that 
there  may  be  no  analytical  solution  for  0(0),  and  numerical  methods  may  need  to  be  applied. 

2  Such  ID  exhaustive  searches  can  also  be  performed  using  tensor-based  approaches  [14]-  [16]. 


Fig.  1.  (a)  Reconstructed  ODF.  (b)  Analytically  defined  ID  space  is  searched,  (c)  All  the 
extrema  and  saddle  points  are  identified,  (d)  ODF  maxima  are  extracted. 


Table  1.  The  associated  Legendre  functions  required  for  the  proposed  algorithm. 


Function 

Expression 
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2.4  One-Dimensional  Exhaustive  Search 


Here  we  detail  the  discretization  scheme  used  to  perform  the  aforementioned  ID 
exhaustive  search  for  the  maxima.  We  exploit  the  closed-form  description  of  the 
curve  (0(0),  0)  provided  by  Eq.  (10)  and  parameterize  the  curve  with  0  6  [0,27i). 
To  achieve  a  constant  spatial  resolution  As  =  AO2  +  sin2  0  A 02,  we  need  a  variable 
step  size  A 0: 


A0 
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1  +  t2(0) 
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which  is  rewritten  as  a  function  of  t(0)  :=  tan  0(0).  For  every  0,  Eq.  (10)  results  in 
one,  two  or  three  real  values  for  t(0),  for  each  of  which  t  (0)  can  be  computed 
simply  by  deriving  Eq.  (10)  with  respect  to  0,  and  substituting  for  0  and  t.  Therefore, 
at  each  step  we  choose  A 0  to  be  the  minimum  of  the  three  (or  fewer)  values  obtained 
from  Eq.  ( 1 1). 

Next,  we  keep  all  the  candidate  points  satisfying  inequalities  (5),  (6),  and  the 
following,  which  is  a  relaxation  of  Eq.  (3), 


I0e  (0,  0)1  <  a. 


(12) 


Fig.  2.  Extracted  maxima  from  synthetic  ODFs  with  fiber  crossing,  in  noise-free  case  (top), 
and  with  SNR=40  (bottom). 


We  found  an  appropriate  value  of  a  =  0.02~0.03  for  the  threshold.  Note  again 
that  inequalities  (5),  (6),  and  ( 12)  can  all  be  computed  analytically  using  equations 
(1),  (2),  (8),  and  Table  1.  The  ODF  maxima  are  then  computed  as  the  mean  directions 
corresponding  to  the  clusters  of  points,  created  by  processing  all  the  candidate  points, 
as  follows:  Each  point  is  added  to  a  previous  cluster  if  its  Euclidean  distance  to  the 
representative  (mean)  point  of  that  cluster  is  minimum  among  all  other  clusters  and  is 
smaller  than  a  threshold  (0.4  was  used  here).  If  no  such  cluster  is  found,  a  new  cluster 
is  created,  and  the  algorithm  goes  on  until  all  the  candidate  points  are  processed. 


3  Results  and  Discussion 

To  validate  our  approach,  we  first  show  results  on  artificial  data.  We  simulated  fiber 
crossing  by  generating  diffusion  images  from  the  sum  of  two  exponentials,  S(u )  = 
(e-“TDlU  +  e~uTD2U)/2,  where  D1  is  a  diagonal  matrix  with  diagonal  entries  (9,  2,  2), 
and  D2  is  D1  rotated  about  the  z-axis  by  a  varying  angle.  CSA-ODFs  were 
reconstructed  in  the  fourth  order  RSSH  basis  from  76  diffusion  directions,  uniformly 
sampled  on  the  sphere.  The  maxima  were  then  extracted  using  the  proposed 
technique,  and  results  are  depicted  in  Fig.  2  (top).  Increasing  the  angular  precision  to 
0.5°  revealed  that  multiple  fiber  orientations  are  resolved  starting  at  the  crossing  angle 
of  37.5°.  Choosing  a  spatial  resolution  of  As  =  0.001,  required  the  evaluation  of  the 
ODF  at  7.7xl04  points,  whereas  a  2D  search  on  the  sphere  with  the  same  resolution 
would  cost  1.6xl07  operations.  When  we  repeated  the  experiment  by  adding  Rician 
noise  with  a  signal-to-noise  ratio  (SNR)  of  40  (Fig.  2,  bottom),  the  minimum  angle 
where  crossing  was  detected  increased  to  48°.  Such  experiments  are  commonly 
employed  to  evaluate  the  robustness  of  the  ODF  reconstruction  algorithm  to  noise. 

We  also  tested  our  method  on  a  popular  public  human  brain  dataset  [17].  CSA- 
ODFs  were  reconstructed  in  the  fourth  order  RSSH  basis  from  200  diffusion  images 
acquired  at  b=3000  s/mm2.  Figure  3  illustrates  the  ODFs  with  their  extracted  maxima 
superimposed  on  the  fractional  anisotropy  (FA)  map,  in  the  region  of  the  centrum 
semiovale,  where  three  major  fiber  bundles  intersect.  To  demonstrate  the  performance 
of  the  proposed  technique,  all  the  maxima  are  shown  here,  including  those 
corresponding  to  slight  variations  in  the  ODF  (for  example  due  to  noise).  Major  ODF 
peaks  corresponding  to  fiber  orientations  may  however  be  selected  by  placing  a 
threshold  on  the  ODF  [3]  or  on  its  curvature  [15], 
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Fig.  3.  Experimental  results  on  human  brain  data,  superimposed  on  the  FA  map. 
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